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We employ a semiclassical picture to study dynamics in a bosonic Josephson junction with vari- 
ous initial conditions. Phase-diffusion of coherent preparations in the Josephson regime is shown to 
depend on the initial relative phase between the two condensates. For initially incoherent conden- 
sates, we find a universal value for the buildup of coherence in the Josephson regime. In addition, 
we contrast two seemingly similar on-separatrix coherent preparations, finding striking differences 
in their convergence to classicality as the number of particles increases. 



I. INTRODUCTION 

Bose-Einstein condensates (BECs) of dilute, weakly- 
interacting gases offer a unique opportunity for exploring 
non-equilibrium many-body dynamics, far beyond small 
perturbations of the ground state. Highly excited states 
are naturally produced in BEC experiments and their 
dynamics can be traced with great precision and control. 
The most interesting possibilities lie in strong correlation 
effects, which imply a significant role of quantum fluctu- 
ations. 

The importance of correlations and fluctuations may 
be enhanced by introducing an optical lattice, that can 
be controlled by tuning its depth. This tight confinement 
decreases the kinetic energy contribution with respect to 
the interactions between atoms. In the tight-binding 
limit, such systems are described by a Bose-Hubbard 
Hamiltonian (BHH), characterized by the hopping fre- 
quency K between adjacent lattice sites, the on-site in- 
teraction strength U, and the total atom number N. The 
strong correlation regime is attained when the character- 
istic coupling parameter u = U N/K exceeds unity, as 
indicated by the quantumphase transition from a super- 
fluid to a Mott-insulator 0, Q • 

The simplest BHH is obtained for two weakly coupled 
condensates (dimer). Its dynamics is readily mapped 
onto a SU(2) spin problem and is closely related to the 
physics of superconductor Josephson junctions [1,01 • To 
the lowest order approximation, it may be described by 
a Gross-Pitaevskii mean-field theory, accurately account- 
ing for Josephson oscillations 041] and macroscopic self 
trapping [9j, observed experimentally in Refs.[l0l.fTl|. as 
well as the equivalents of the ac and dc Josephson effect 
observed in [isj . 

Both Josephson oscillations and macroscopic self trap- 
ping rely on coherent (Gaussian) preparations, with 
different initial population imbalance. The mean-field 
premise is that such states remain Gaussian throughout 
their evolution so that the relative phase if between the 
two condensates remains defined. However, interactions 



between atoms lead to the collapse and revival of the rel- 
ative phase in a process known as phase diffusion 
(the appropriate term is in fact phase spreading). Phase 
diffusion has been observed with astounding precision in 
an optical lattice in Refs. 17|, in a double-BEC system in 
Refs.[18 20], and in a ID spinor BEC in Ref. [21]. Typ- 
ically, the condensates are coherently prepared, held for 
a varying duration ('hold time') in which phase-diffusion 
takes place, and are then released and allowed to inter- 
fere, thus measuring the relative coherence through the 
many-realization fringe visibility. In order to establish 
this quantity, the experiment is repeated many times for 
each hold period. 

Phase-diffusion experiments focus on the initial prepa- 
ration of a zero relative phase and its dispersion when no 
coupling is present between the condensates. However, in 
the presence of weak coupling during the hold time, the 
dynamics of phase diffusion is richer. It becomes sen- 
sitive to initial value of (p and the loss of coherence is 
most rapid for <y9 = tt (2342^ . Here, we expand on a re- 
cent letter ^2^], showing that this quantum effect can be 
described to excellent accuracy by means of a semiclas- 
sical phase-space picture. Furthermore, exploiting the 
simplicity of the dimer phase space, we derive analytical 
expressions based on the classical phase-space propaga- 
tion [H. 

Phase-space methods [265 have been extensively ap- 
plied for the numerical simulation of quantum and ther- 
mal fluctuation effects in BECs |27h39| . Such meth- 
ods utilize the semiclassical propagation of phase-space 
distributions with quantum fluctuations emulated via 
stochastic noise terms, and using a cloud of initial condi- 
tions that reflects the uncertainty of the initial quantum 
wave-packet. One particular example is the truncated 
Wigner approximation [H, [sil, IH, |3g| where higher order 
derivatives in the equation of motion for the Wigner dis- 
tribution function are neglected, thus amounting to the 
propagation of an ensemble using the Gross-Pitaevskii 
equations. 

Due to the relative simplicity of the classical phase- 
space of the two-site BHH, it is possible to carry out its 
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semiclassical quantization semi- analytically [H, [H, g^- 
H^l and acquire great insight on the ens uing dynamics of 
the corresponding Wigner distribution [4J, |43| • In this 
work, we consider the Jopheson-regime dynamics of four 
different preparations, interpreting the results in terms 
of the semiclassical phase-space structure and its impli- 
cations on the expansion of each of these initial states in 
terms of the semiclassical eigenstates. We first explore 
the phase-sensitivity of phase-diffusion in the Josephson 
regime [H, . Then we study the buildup of coherence 
between two initially separated condensates [4l| , which is 
somewhat related to the phase-coherence oscillations ob- 
served in the sudden transition from the Mott insulator 
to the superfluid regime in optical lattices jisl - liTj . Then 
we compare two coherent preparations in the separatrix 
region of the classical phase space, finding substantial 
differences in their dynamics due to their different par- 
ticipation numbers. 

The narrative of this work is as follows ja] : The two- 
site Bose-Hubbard model and its classical phase-space are 
presented in Section II. The semiclassical WKB quantiza- 
tion is carried out in Section III. In Section IV, we define 
the initial state preparations of interest and explain their 
phase-space representation. This is used in Section V to 
evaluate their expansion in the energy basis (the local 
density of states), which is the key for studying the dy- 
namics in Sections VI to VIII. Both the time-averaged 
dynamics and the fiuctuations of the Bloch vector are 
analyzed. In particular we observe in Section VIII that 
the fiuctuations obeys a remarkable semiclassical scaling 
relation. Conclusions are given in Section X. 



conserving the spin = i(j + 1) with j = N/2. The 
BHH Hamiltonian thus has a spherical phase-space struc- 
ture. In the absence of interaction ([/ = 0) it describes 
simple spin precession with frequency i7 = (K,0,£). 

The quantum evolution of the spin Hamiltonian (j4]) is 
given by the unitary operator exp{~i'Ht). Its classical 
limit is obtained by treating {Jx, Jy, Jz) as c-numbers, 
whose Poisson Brackets (PB) correspond to the SU(2) 
Lie algebra. The classical equations of motion for any 
dynamical variable A are derived from A = — [H, A]pb. 
Conservation of allows for one constraint on the non- 
canonical set of variables (Jx, Jy, Jz)- 

Due to the spherical phase-space geometry, it is natural 
to use the non-canonical conjugate variables (tp, 0), 



Jz 

J X 



[J2]l/2cos(0) 

sin(6/) cos(<^) 



(5) 
(6) 



where [ J^]^/^ is a constant of the motion. In the quantum 
mechanical Wigner picture treatment (see later sections) 
this constant is identified as [(j+l)j]^^^, based on the 
derivation of section 2.3 of Ref. [i^. For large j we use 
w N/2. With these new variables the Hamilto- 
nian takes the form 
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where the scaled parameters are 
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II. THE TWO-SITE BOSE-HUBBARD MODEL 

We consider the Bose-Hubbard Hamiltonian for N 
bosons in a two-site system. 
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where K is the hopping amplitude, U is the interaction, 
and £ — £2 — £i is the bias in the on-site potentials. We 
use boldface fonts to mark dynamical variables that are 
important for the semiclassical analysis, and use regu- 
lar fonts for their values. The total number of particles 
ni -\- n2 — N is conserved, hence the dimension of the 
pertinent Hilbert space is A/" = iV + 1. Defining 
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J+ = a\a2 



(2) 
(3) 



and eliminating insignificant c- number terms, we can re- 
write Eq. ([1} as a spin Hamiltonian, 



Ur-£J,-KJx, 



The structure of the underlying classical phase space is 
determined by the dimensionless parameters u and e. 
For strong interactions u > 1 the phase space includes a 
figure-eight shaped separatrix (Fig.[T]), provided |e| < Sc, 
where Um, El 



(9) 



This separatrix splits phase space into three inte- 
grable regions: a 'sea' of Rabi-like trajectories and two 
interaction-dominated nonlinear 'islands'. For zero bias 
the separatrix is a symmetric 8 shaped figure that en- 
closes the two islands, and the relevant energies are [i^ 



= 

E^ = 



-{l/2)NK 
+ {l/2)NK 



E+= {l/'\)[u+{l/u)\NK = top energy. 



— ground energy, (10) 
= separatrix, (11) 
(12) 



Looking at the phase-space structure as a function of u, 
the two islands emerge once u > 1. For w > 2, they en- 
compass the North and the South poles, yielding macro- 
scopic self-trapping [9] , and for very large u ^ 1 they 
cover most of the Northern and the Southern hemi- 
spheres respectively. Note that for u ^ 1, the expres- 
sion E+ K, {N/2)'^U, reflects the cost of localizing all the 
(4) particles in one site, compared to equally-populated sites. 



3 



As an alternative to Eq. ([T]) , it is possible to employ the 
relative number-phase representation, using the canoni- 
cally conjugate variables (tp, n), with the Hamiltonian, 



n = Un^ - £n - Ky/iN/2y - cos{ip). (13) 

In the \n\ ^ {N/2) region of phase space, one obtains the 
Josephson Hamiltonian, which is essentially the Hamil- 
tonian of a pendulum 



= Ec{n-n,f - Ejcos{cp) (14) 



with Ec — U and Ej = KN/ 2 while is linearly related 
to £. The Josephson Hamiltonian ignores the spherical 
geometry, and regards phase space as cylindrical. Ac- 
cordingly it captures much of the physics if the motion 
is restricted to the small-imbalance slice of phase-space, 
as in the case of equal initial populations and strong 
interactions. However, for other regimes it is an over- 
simplification because it does not correctly capture the 
global topology. 



III. THE WKB QUANTIZATION 

We begin by stipulating the procedure for the semiclas- 
sical quantization of the spin spherical phase-space. In 
the {(fi, n) representation, the area element is dO. = dipdn, 
so that the total phase-space area is 27rN with Planck 
cell h = 27r. Alternatively, using the {ip, 9) coordinates, 
the area element is dil = dipd cos 6, so that the total area 
is 47r. Consequently 



Planck cell area in steradians 



An 

77' 



(15) 



Within the framework of the semiclassical picture, eigen- 
states are associated with stripes of area h that are 
stretched along contour lines H{ip,d) — E oi the classi- 
cal Hamiltonian 7i. The associated WKB quantization 
condition is 



A{E,) 



0,1,2, 



(16) 



where A{E) is defined as the phase-space area enclosed 
by an E contour in steradians: 



A{E) = Q{E-n{^,9)) dn, 



(17) 



while the area of phase space in Planck units is A{E)/h. 
The frequency of oscillations at energy E is 



ujiE) 



dE 
dv 



A'{E) 



(18) 



Consider the semiclassical quantization of the Hamil- 
tonian ([7]). We distinguish three regimes of interaction 
strength. Assuming e = 0, these are 

Rabi regime: u < 1 (no islands) (19) 
Josephson regime: 1 < u < N"^ (see Fig. [T]) (20) 
Fock regime: u > N"^ (empty sea) (21) 



In the Fock regime, the area of the sea becomes less than 
a single Planck cell, and therefore cannot support any 
eigenstates. Our interest throughout this paper is mainly 
in the Josephson regime where neither the K term nor 
the U term can be regarded as a small perturbation in 
the Hamiltonian. This regime, characteristic of current 
atom-interferometry experiments, is where semiclassical 
methods become most valuable. 

In the WKB framework the spacing between energies 
equals the characteristic classical frequency at this en- 
ergy. If the interaction is zero {u — 0), the energy lev- 
els are equally spaced and there is only one frequency, 
namely the Rabi frequency 



UJK = K 



(22) 



For small interaction (0 < m < 1) the frequencies around 
the two stable fixed points are slightly modified to the 
plasma frequencies: 



= u{E±) = y/{K±NU)K. 



(23) 



If the interaction is strong enough {u ^ 1), the oscillation 
frequency near the minimum point can be approximated 
as 



ujj = u}{E^) w VNUK = ^UJK, (24) 

while at the top of the islands we have: 

a;+ = a;(£'+) w NU ^ u ojk- (25) 

The associated approximations for the phase-space area 
in the three energy regions (bottom of the sea, separatrix, 
and top of the islands) are respectively 



\AiE) = \AiEM 
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(26) 
(27) 
(28) 



In the last expression the total area of the two islands 
{E^~E) /U should be divided by two if one wants to 
obtain the area of a single island. A few words are in 
order regarding the derivation of Eq. (|27|) . In the vicin- 
ity of the unstable fixed point the contour lines of the 
Hamiltonian are C/n^ — {NK/A)ip^ = E—E^, where for 
convenience the origin is shifted {ip=n i— 0). Defining 
= 4{E~E^) /NK the area of the region above the sep- 
aratrix is 2{NK/my^[A{a) - A{0)] where A{a) is the 
integral over y^(p^ + . The result of the integration 
is A{a) — A{0) = log(5/a), with some ambiguity with 
regard to 6 ~ 1 which is determined by the outer limits 
of the integral where the hyperbolic approximation is no 
longer valid. 

Away from the separatrix, the WKB quantization 
recipe implies that the local level spacing at energy E 
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is u}{E) given by Eq. ([18]). In particular, the low energy 
levels have spacing uij, while the high energy levels are 
doubly-degenerate with spacing uj+. In the vicinity of 
the separatrix we get 





NK 




-log 






TT 


E ~ Ex 





Using the WKB quantization condition, we find that the 
level spacing at the vicinity of the separatrix {E ^ E^) 
is finite and given by the expression 

Lo^ = [\og{N/^)Y^ u:j. (30) 

Using an iterative procedure one finds that at the same 
level of approximation the near-separatrix energy levels 
are 



E^^E^ + 


-log 








TT 


v—v^ 





{v-v^ wj, (31) 



where = A{E^ / h. Fig.[T]dcmonstrates the accuracy of 
the WKB quantization, and of the above approximations. 

IV. THE INITIAL PREPARATION AND ITS 
PHASE-SPACE REPRESENTATION - THE 
WIGNER FUNCTION 

Our approach for investigating the dynamics of various 
initial preparations relies on the Wigner-function formal- 
ism for spin variables, developed in Refs. [isl. |4^. Each 
initial preparation is described as a Wigner distribution 
function over the spherical phase space. The dynam- 
ics is deduced from expanding the initial state in terms 
of the semiclassical eigenstates described in Sec III. In 
this section we specify the Wigner distribution for the 
preparations under study whereas the following section 
presents the eigenstate expansion of each of these four 
initial wavepackets, evaluated semiclassically. 

To recap the phase-space approach to spin [H, [l^l , the 
Hilbert space of the BHH has the dimension M = 2j-|-l, 
and the associated space of operators has the dimension- 
ality A/"^. According to the Stratonovich- Wigner- Weyl 
correspondence (SWWC) [s^l, any observable A in this 
space, as well as the probability matrix of a spin'^-'^ entity, 
can be represented by a real sphere'^-'^ function Ay^iyi). 
The sphere'^^-') is spanned by the Y^^i^) functions with 
£ <2j, and the practical details regarding this formal- 
ism can be found in Refs. [H, [3. The SWWC allows 
to do exact quantum calculation in a classical-like man- 
ner. A few examples for Wigner functions pertinent to 
this work, are displayed in Fig. [5] Expectation values are 
calculated as in classical statistical mechanics: 

—p„{n)A„{n) (32) 

In particular the Wigner- Weyl representation of the 
identity operator is 1, and that of Jx is as expected 



sm{d) cos{ip) 0]. We adopt the convention 
that is normalized with respect to the measure dfl/h, 
allowing to handle on equal footing a cylindrical phase 
space upon the re-identification dQ — dipdn and h ^ 2tt. 

Within this phase-space representation, the Fock 
states \n) are represented by stripes along constant 9 con- 
tours (see e.g. left panel of Fig. [5]). The \n=N) state (all 
particles in one site) is a Gaussian-like wave packet con- 
centrated around the NorthPole. From this state, we can 
obtain a family of spin coherent states (SCS) \0,(p) via 
rotation. 

In what follows, we explore the dynamics of the follow- 
ing experimentally-accessible preparations (see Fig. [2]), 
the first being a Fock state, whereas the last three are 
spin coherent states: 

• TwinFock preparation: The n=0 Fock preparation. 
Exactly half of the particles are in each side of the 
double well. The Wigner function is concentrated 
along the equator 6 — tt/2. 

• Zero preparation: Coherent {0=^/2, ip=0) prepa- 
ration, located entirely in the (linear) sea region. 
Both sites are equally populated with definite rel- 
ative phase. The minimal wave-packet is centered 
at (n = 0, iy9 = 0). 

• Pi preparation: Coherent {6—T:/2,ip—Tr) on- 
separatrix preparation. The sites are equally pop- 
ulated with TT relative phase. The minimal wave- 
packet is centered at (n = 0, = tt). 

• Edge preparation: Coherent 95 7^ tt on-separatrix 
preparation. The minimal wave-packet is centered 
on the separatrix but away from the saddle point 
on the ip—O side. 

The Wigner function of an SCS resembles a minimal 
Gaussian wave-packet, and it should satisfy 

j p^{6M^^ J[p„ie,^)rf^l. (33) 

This requirement helps to determine the phase space 
spread without the need to use the lengthy algebra of 
Refs. [H, Q . For the Fock coherent state \n=N) , that 
is centered at the NorthPole (9 = 0), one obtains: 

p^^\9,^)^2e-^'' . (34) 

For the coherent states centered around the Equator, it is 
more convenient to use (ip, n) coordinates, e.g. the f — 
coherent state is well approximated as, 

p\^\n,(p) w i-e^^^ip- , (35) 
ab 

with a = 1/VW and b ^ ^^71/2. Shifted versions of 
these expressions describe the ip = tt and the "Edge" 
preparations. The Wigner function of a Fock state 
"0 = |n) is semiclassically approximated as 

Pw\n,ip) w (5(n-n), (36) 
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whereas the Wigner function of an eigenstate is semiclas- 
sically approximated by a micro-canonical distribution: 



(37) 



In general none of the above listed initial states is an 
eigenstate of the BHH, but rather a superposition of BHH 
eigenstates. Consequently their Wigner function deforms 
over time (see e.g. Fig. [3]) and the expectation values of 
observables become time dependent (see e.g. Fig. |4|), as 
discussed in later sections. 



V. THE INITIAL PREPARATION AND ITS 
EIGENSTATE EXPANSION - LOCAL DENSITY 
OF STATES 

Having set the stage by defining the phase-space rep- 
resentation of eigenstates in Sect. Ill and of the initial 
conditions in Sect. IV, the evolution of any initial prepa- 
ration is uniquely determined by its eigenstate expansion. 
Thus, in order to analyze the dynamics, we now evaluate 
the local density of states (LDOS) with respect to the 
preparation -0 of the system: 



V{E,) = |(ii;,|V^)p =trace(pMpW) 



= P 



(38) 



If ip is mirror symmetric the above expression should be 
multiplied either by or by 2 in the case of odd/even 
eigenstates. In Fig. [5l we plot the LDOS associated with 
Pi, Edge, Zero and TwinFock preparations. The applica- 
bility of semiclassical methods to calculate the LDOS has 
been numerically demonstrated for the case of a three site 
(trimer) model in [s^. By contrast, the simpler two site 
(dimer) model under consideration, offers an opportunity 
for deriving exact expressions by substituting Eqs. (j35p - 
(|37)) into Eq. ((38|) and evaluating the integrals under the 
appropriate approximations. For clarity, we summarize 
below the main results of this analytic evaluation, with 
the details given in App. [D] 

Consider first the Pi and Edge preparations. The 
Wigner distributions of both lie on the separatrix and 
are hence concentrated around the energy E^. Yet, their 
line shapes are strikingly different: For an Edge prepa- 
ration we have, 



Y>{E) oc ^ exp 



1 



E-E^ 



(39) 



featuring a dip dX E ^ E^. The calculation in the Pi case 
leads to 



P(S) (X '^^^^ Bessel 



E-E^ 
NU 



(40) 



(the exact expression can be found in App. [D]). This 
latter line shape features a peak at E ^ E^, because the 



Bessel function compensates the logarithmic suppression 
by uj{E). Consequently, as seen below, the fluctuations 
associated with on-separatrix motion differ dramatically, 
depending on where the SCS wave-packet is launched. 

Analytic results are obtained in App. |D1 also for the 
Zero and TwinFock preparations. In the latter case the 
result is: 



P{E) cx 



NK 



2E 

'nk 



-1/2 



(41) 



As shown in Fig. [SI the above semiclassical expressions 
agree well with the exact numerical results. 

The LDOS determines the spectral content of the dy- 
namics. Consider first the characteristic oscillation fre- 
quency of dynamical variables. Away from the separa- 
trix, it is given by the classical estimate u)^^^ ~ '^{E)- 
For example, for a Zero preparation, we have Wo^^ = wj. 
While the classical frequency vanishes on the separtrix, 
we still obtain a finite result for near-separatrix prepa- 
rations (Pi, Edge) because the wave-packet has a finite 
width, and because uj^ provides a lower bound on lo^^^. In 
any case the result becomes h dependent. To be specific, 
consider the Pi and the Edge preparations separately. In 
both cases the width of the wave-packet is An = N/2. 
In case of the Pi preparation it occupies an energy range 
A£' = [/An^ (X N , while in the Edge preparation case 
IS.E = t;„An cx N^^^, where Vn ~ ujj is identified as the 
velocity of the phase space points in this region. Using 
Eq. ^ we find 



1 ( ^ 
log — 

u 



(42) 



where the additional factor of "2" applies to the Edge 
preparation: This factor is due to the different depen- 
dence of AE on N in the two respective cases. On top 
we might have an additional factor of 2 due to mirror 
symmetry (see discussion of the dynamics at the end 
of the next section). Note that as u/N exceeds unity, 
the distinction between Pi, Zero, and Edge preparations 
blurs because the wave-packets becomes wider than the 
width of the separatrix region, until at the Fock regime, 
all three consist of the same island levels. Thus in this 
limit, we get for Wosc essentially the same result as in the 
(negligable K) Fock regime: the width of the wave-packet 
is An = \/nJ2, and the dispersion relation lij{E) = 2[/n 
gives the standard separated-condensates phase-diffusion 
frequency [13, [11, 



/ M \ 1/2 I 



(43) 



It should be clear that, both, in Eq. (g^]) and in Eq. P5)) 
u/N should be accompanied with a numerical prefactor 
that should be adjusted because the notion of "width" is 
somewhat ill defined and in general depends on the pre- 
cise details of the numerical procedure, which we discuss 
in the next section. 
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For the analysis of the temporal fluctuations it is cru- 
cial to determine the number of eigenstates that partici- 
pate in the wave-packet superposition. This is given by 
the LDOS participation number, 



M 



(44) 



See Fig. [6] for numerical results. The participation num- 
ber can be roughly estimated as M = AE/uta^^, where 
AE is the energy width of the wave-packet and lu^^^ is 
the mean level spacing. In the case of a Pi preparation 









M ; 




h(v)] 









while for an Edge preparation we find 



M 



log 











N. 



(45) 



(46) 



Note that M/viV is a function of the semiclassical ra- 
tio (N/uy^^ between the energy width of the wave- 
packet and the width of the separatrix region. The ex- 
pected scaling is confirmed by the numerical results of 
Fig.m The above approximations for M assume u/N < 1 
and are useful for the purpose of rough estimates. 

In the next sections we analyze the temporal fiuctua- 
tions of some observables. The associated Fourier power 
spectrum (Fig. [Sj right panels) is related to the LDOS 
content of the wave-packet superposition. Both the char- 
acteristic frequency (Fig. [7]) and the spectral spread of 
the frequencies can be estimated from the Wo^c and the 
M that are implied by the above LDOS analysis. 



VI. DYNAMICS (I) - THE TIME EVOLUTION 
OF THE BLOCH VECTOR 



After describing the semiclassical phase-space picture 
and using it to determine the spectral structure of the ini- 
tial preparations, we now turn to the ensuing dynamics. 
At present, most experiments on the bosonic Josephson 
system, measure predominantly quantities related to the 
one-body reduced probability matrix, defined via the ex- 
pectation values Si = {2/N){Ji) as. 



1 ,.t~ ^ 



(47) 



where S = {Sx,Sy, Sz) is the Bloch vector and & is com- 
posed of Pauli matrices. The population imbalance and 
relative phase between the two sites, determined by di- 
rect imaging and the position of interference fringes , 
are given respectively by, 

OccupationDiff = N S^, (48) 
RelativePhase = arctan(S'3;, 5'j,), (49) 



whereas single-particle purity is reflected by the mea- 
sures, 

OneBodyPurity ^ (1/2) [l+S^+S^+S^,] , (50) 
Fringe Visibility = [S^ + S^^^^ . (51) 

In particular the Fringe Visibilty, given by the transverse 
component of the Bloch vector, corresponds to the visi- 
bility of fringes, averaged over many realizations [l9l[20|. 
Loosely speaking it reflects the phase-uncertainty of the 
state. Coherent states have maximum OneBodyPurity. 
Of these, equal-population coherent states have maximal 
Fringe Visibility with the smallest phase-variance. Start- 
ing from a coherent preparation, single-particle purity 
can be diminished over time due to nonlinear effects (see 
below) or due to interaction with environmental degrees 
of freedom (decoherence). By contrast, Fock states carry 
no phase information, but interactions may lead to their 
dynamical phase-locking over time and to the buildup of 
Fringe Visibility (see below). 

We carry out numerically-exact quantum simulations 
where the state \tp) is propagated according to \ipit)) = 
exTp{—it'H)\^) , where the Hamiltonian is given by Eq.Q 
which is equivalent to Eq.(IT]). The evolved state after 
time t can be visualized using its Wigner function. For 
example, the time evolution of the initial TwinFock state 
is illustrated in Fig. [31 with the Wigner function of an 
evolved state shown in Fig. [Sj^a). Comparison is made 
in Fig. [3Kb) to the Liouville propagation for the same 
duration, of the corresponding cloud of points, according 
to the classical equations. 



Sy 



-(1 + uSx)Sz 



(52) 
(53) 
(54) 



where time has been rescaled (t := Kt). Good quantuni- 
to-classical correspondence is observed for short-time 
simulation (e.g. see Fig. [J]). In Fig. [SJc) we plot the 
resulting occupation statistics, which can be regarded as 
a projection of the phase space distribution (classical, 
dash-dotted line) or Wigner function (quantum, solid 
line), namely Pt(n) = \{n\^p{t))\^ = trace(p(")p('^(*»). 
Our main interest is in the Fringe Visibility, and hence in 

S. it) = ^ ( J.) = '.^ (sin(0) cos(v)) . (55) 

N J 

The prefactor in the last equality is implied by Eq. ([5|), 
and cannot be neglected if the number of particles is 
small. In Fig. [4] we plot Sx{t) for the four prepara- 
tions defined in Sec. IV, comparing semiclassical results 
(dash-dotted lines) to the numerical full quantum cal- 
culation (solid lines). For all equatorial preparations 
(Zero, Pi, TwinFock), Sx is in fact the Fringe Visibility, be- 
cause Sy{t) — identically throughout the evolution. As 
a general observation, the semiclassical simulation cap- 
tures well the short-time transient evolution and the long 
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time average Sx- For example, for the initial Twin- 
Fock preparation (Fig. HJi), it reproduces the univer- 
sal Josephson-regime Fringe Visibility of ~ 1/3, resulting 
from the dynamical smearing of the Wigner distribution 
function throughout the linear sea region of phase-space 

In what follows we would like to determine the long 
time average Sx and the power spectrum C{uj) = 
FT[f{t)] of the temporal fluctuations f{t) = Sx{t) - S!^. 
(FT denotes Fourier Transform). Characteristic power- 
spectra for the pertinent preparations are shown in the 
right panels of Fig. [S] We characterize the fluctuations 
by their typical frequency w^ac, by their spectral support 
(discrete or continuous- like) , and by their RMS value: 



RMS[Sx] = [f{tY 



1/2 _ 



C{uj)duj 



1/2 



(56) 



The dependence of and Sx, and RMS [5*1,] on the 
dimensionless parameters (u, N) is illustrated in Figs.[7]- 
m It should be realized that the observed frequency of 
the Sx{t) oscillations is in fact 2uj^^^ due to the mirror 
symmetry of the observable. 



VII. DYNAMICS (II) - FRINGE VISIBILITY IN 
THE JOSEPHSON REGIME 

In the Fock regime, the Fringe Visibility of an ini- 
tial coherent preparation decays to zero: the ini- 
tial Gaussian-like distribution located at (6'=7r/2,(^) is 
stretched along the equator, leading to increased relative- 
phase uncertainty with fixed population-imbalance. This 
phase spreading pro cess is known in the literature as 
'phase diffusion' [l4| - [l7l . [2l|. By contrast, a TwinFock 
preparation is nearly an eigenstate of the BHH in this 
regime, and its zero Fringe Visibility remains vanishingly 
small from the beginning. 

In the Josephson regime the dynamics of the single- 
particle coherence is more intricate. In this section we 
discuss the evolution of S{f) within the framework of 
the semiclassical approximation. As a first step, we dis- 
regard fluctuations and recurrences and address only the 
time- averaged dynamics. In particular, we determine the 
long-time average of the Fringe Visibility, which f or the 
TwinFock, Pi, and Zero preparations is given by Sx{t), 
as noted after Eq. (l55t . 

Numerical results for the long-time average and for the 
RMS of the fluctuations as a function of u/N are pre- 
sented in Fig. [5] and further discussed below. The main 
observations regrading the dynamics in the Josephson 
regime are: 

• The TwinFock preparation of fully-separated con- 
densates devel ops p hase-locking at (p^O, with 
FringeVisibility Sx{t) « 1/3 [U. 

• Starting from a SCS preparation, phase-diffusion 
becomes phase sensitive. The Zero preparation is 



phase-locked, while the coherence of the Pi prepa- 
ration is partially lost [s^] , exhibiting huge 
fluctuations. 

• The Edge preparation exhibits distinct behavior, 
that neither resembles the Zero nor th e Pi prepa- 
rations, involving sign reversal of Sx{t). 

In the remaining part of this section we quantify these 
observations by flnding the long-time average Sx (t) based 
on simple phase-space considerations. In the Joseph- 
son regime, the value of Sx{t) for a coherent prepa- 
ration should be determined by the ratio between its 
An w N/2 width and the width of the separatrix re- 
gion An « y/NK/U, i.e. 



the semiclassical ratio = {u/Ny^'^ 



(57) 



This ratio determines the long-time phase-space distribu- 
tion of the evolving semiclassical cloud: In the case of a 
TwinFock preparation this cloud fills the entire sea re- 
gion; in the " Zero" case it is confined to an ellipse within 
the sea region; and in the "Pi" case it stretches along 
the separatrix and therefore resembles a micro-canonical 
distribution. The projected phase-distribution P{(p) is 
determined accordingly. Disregarding a global normal- 
ization factor we get 



p(^) 
p(^) 



exp[-(^V(4w/^)] [Zero] , (58) 
[(u/iV)-|-cos2(^/2)]-V2 [Pi] , (59) 

cos(.^/2) [TwinFock] . (60) 



A few words are in order regarding the derivation 
of the above expressions. The Zero case prepa- 
ration is represented by the Gaussian of Eq. ([55)) 
whose major axis is 

An « (7V/2)i/2. This Gaussian 
evolves along the contour lines of the Hamiltonian 
H(n, (p) = Un^ + {NK/2) cos{p). After sufficiently long 
time the evolving distribution still has the same An, 
but because of the spreading its other major axis, as 
determined by the equation UAri^ = {NK/A)A(p'^, be- 
comes A.^ « (2f//X)i/2 leading to Eq.dSH]). The Twin- 
Fock preparation is represented by Eg. ([55)1 . After suffi- 
ciently long time the evolving distribution fills the whole 
sea H(0, (f) < -Ex- The equation that describes this filled 
sea can be written as n < n^{ip), where 



' NK 

"x(<y5) = \l (1 + C0S((^)) . 



(61) 



The projection of area under n < n^{(p) is simply 
P{ip) oc n^{ip), leading to Ea. (l5D|) . The Pi case prepa- 
ration is represented by the Gaussian that is located on 
the separatrix. After sufficiently long time the evolv- 
ing distribution is stretched along n ~ n^ip), and looks 
like d{'H{6, ip) — E^). If we neglected the finite width of 
this distribution we would obtain P{ip) oc l/n^[(p), which 
is divergent at ip t:. But if we take into account the 
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An \/N/2 width of the preparation, which is effec- 
tively like adding {N/2) under the square root of Eq. ifGl]) . 
then we get Eq. lfSQ]) . 

As implied by the Wigner-Weyl picture one can get 
from P(</9) the long time average Sx through the integral 



cos((^) P{(p)d(p 



(62) 



The calculation is straightforward leading to 

exp[-{u/N)] [Zero], (63) 



Sx^ -1 -4/ log [^{u/N) 
1/3 



[Pi], (64) 
[TwinFock]. (65) 



These expressions agree with the numerics of Fig. [51 and 
confirm the predicted scaling with u/N. 



VIII. DYNAMICS (III) - LONG TIME 
FLUCTUATIONS 

To complete our phase-space characterization of the 
Bloch vector dynamics, we need to address the long-time 
fluctuations f{t) from the average value Sx- As demon- 
strated in Fig. m the evaluation of Sx could be done to 
excellent accuracy based on the semiclassical propagation 
of phase-space distributions according to purely classical 
equations of motion. This corresponds to the truncated 
Wigner approach of quantum optics [H, [U 113, [SBj , re- 
taining only the leading Liouville term in the equation of 
motion for the Wigner distribution. Obviously one can 
not guarantee that the remaining Moyal-bracket terms, 
which are initially 0{1/N), will remain small through- 
out the evolution. This is the source of the fluctuations 
observed in Fig. S) While their characterization goes be- 
yond the lowest-order truncated Wigner semiclassics, we 
still can estimate them based on the phase-space LDOS 
expansion, as described below. 

The two-site BHH has essentially one degree of free- 
dom since both the energy and the number of the bosons 
is conserved. Therefore, away from the separatrix, level 
spacing is determined by the classical frequency ijj{E) 
with small h dependent corrections. This should be 
contrasted with higher dimensions d > 1, for which the 
level spacing a h'^~^ is highly non-classical. For d — 1, 
the only region where h determines the level spacing 
oc |log(ft.)|~^ is in the vicinity of the separatrix as im- 
plied by Eq. ([50]). 

The Heisenberg time is defined as the inverse of the 
mean spacing of the participating levels. For a d = 1 sys- 
tem and away from the separatrix, the Heisenberg time 
is merely the period of classical oscillations. Thus in the 
case of a Zero preparation we have uj^,^ — ujj (the ob- 
served frequency is doubled due to Mirror symmetry). 
Close to the separatrix, uj^^^ becomes h dependent as in 
Eas. (|42j|43p . This prediction is confirmed by the numer- 
ics (see Fig. [7]), including the non-symmetry related fac- 
tor 2 that distinguishes the Edge from the Pi preparation. 



Assume the system is prepared in some state V': 6.g. 
a Gaussian-like SCS. We define M as in Eq. (gl]), im- 
plying that V is roughly a superposition of M energy 
states. If the energy levels are equally spaced the motion 
is strictly periodic. Otherwise it is quasi-periodic. Our 
aim is to trace the non-classical behavior in the RMS of 
the temporal fluctuations of an observable A, say of the 
Fringe Visibility as defined in Eq. (j56p . 

Before proceeding, it should be made clear that the 
RMS of the fluctuations of any observable A in a classi- 
cal simulation (i.e. classical propagation of a single tra- 
jectory) is non-zero and characterized by its power spec- 
trum Cci(aj). However, the RMS of the fluctuations in the 
semiclassical evolution (i.e. classical, leading-order prop- 
agation of a cloud of trajectories emulating the Wigner 
function) goes to zero due to the ergodic-like spreading 
of the wave-packet. In contrast, the RMS of the fluctua- 
tions in the quantum evolution (corresponding to the full 
propagation to all orders, of the Wigner distribution) de- 
pends on M. This dependence on M can be figured out 
by expanding the expectation value in the energy basis 



A e''(^- 



(66) 



where ij^u — {E^\il)). The time average of this expectation 

P(i;,)ofEq. 



value is {A)^ ^^ijPvA^^i, where 
This average has a well-defined ^.-independent classical 
limit. But if we first square, and then take the time 
average we get 



(67) 



The matrix elements can be evaluated semiclassically us- 
ing the well know relation = C^-i{E^—E^)/{2tiqY 
where g is the mean level spacing (see Eq.(6) of Ref.^SlTj 
and references therein). For presentation purpose it is 
convenient to visualize C^iiuj) as having a rectangular- 
like lineshape of width w^,, such that its total area is 
C(0) X LOc- It is also convenient to define the dimension- 
less bandwidth b as the spectral width of Cci(aj) divided 
by the mean level spacing, namely h = quJc- 

Using the semiclassical estimate for the matrix ele- 
ments, we consider the outcome of Ea. (j67p . in two lim- 
iting cases. If the energy spread of the wave-packet is 
smaller than the spectral bandwidth, we can factor out 
C^i(0)/(27rf3), and carry out the summation ^p^P^l — 1, 
leading to 



RMS [{A)^] 



Cci{uj)duj 



-.1/2 



(68) 



For integrable one-dimensional systems 6^1 reflects 
that only nearby levels are coupled. A semiclassically 
large bandwidth b oc h^~'^ is typical for chaotic systems, 
which is not the case under consideration. Therefore 
we turn to the other possibility, in which the energy 
spread of the wave-packet is large compared with the 
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spectral bandwidth. In such case we can make in Eq. ((67)) 
the replacement p^, i— >■ 1/M, and consequently the sum 
J2 I^i'.mP equals M times the area of C^j{uj), leading to 



RMS[(A),] = 



— / C4uj)duj 



1/2 



(69) 



This is the same as the classical result but suppressed by 
factor 1/^/M. Note again that the semiclassical result is 
always zero, and corresponds formally to M = c». 

Consider now the RMS of {t) . For TwinFock prepa- 
ration it follows from Eq. (|4T|) that M (x N and hence 
Eq.(|Ml) implies l/N^^^ suppression of the RMS. For co- 
herent preparations Cci(aj) becomes N dependent too, 
and consequently from the discussion after Eq. (|46p it 
follows that the RMS is a function of the semiclassical ra- 
tio Eq. ([57]) . and multiplied by l/A^^/"* suppression factor 
that spoils the semi-classical scaling. This is confirmed by 
the numerics (Fig. [5]). If the dynamics is very close to the 
separatrix the classical fluctuations are 0(1) and there- 
fore the quantum result is RMS [S.j:{t)] « 1/Vm. The 
implication for the on-separatrix coherent preparations. 
Pi versus Edge, is striking: Substitution of Eas. (|45ll46)) 
into Eq. (|69| leads to 



RMS [S^{t)] 



TV- for Edge 
(log(iV))-i/2 for Pi 



(70) 



Thus, convergence to classicality is far more rapid for the 
Edge-preparation than it is for the Pi-preparation, even 
though both lie on the separatrix. With the Pi prepara- 
tion, even if TV is very large (small "/i"), quantum fluc- 
tuations still remain pronounced. In fact, from Eq. ps)) 
it follows that the fluctuations in the Pi case are mainly 
sensitive to the strength u of the interaction. 

Finally, we mention that the analysis of fluctuations 
above is somewhat related to the discussion of thermal- 
ization in Ref. |53i] . and we would like to further connect 
it with the observation of collapses and quantum revivals 
as discussed e.g. in Ref. . Relating to the LDOS, as 
defined in Section V, we note that the collapse time is the 
semiclassical time which is determined by the the width 
of the classical envelope, while the revival time is related 
to the spacing between the spectral lines. The latter can 
be calculated using the formula 



(71) 



with the WKB estimate for i?^ in Section III, leading in 
the separtatrix region to ^ 2t:/uj^. 



IX. CONCLUSIONS 



WKB quantization. Thus, closed semiclassical results 
are obtained for the local density of states of the vari- 
ous initial preparations, providing useful insights for the 
associated quantum evolution. 

Within the framework of mean- field theory (MET) , the 
dynamics is obtained by evolving a single point in phase 
space, using the Gross-Pitaevskii (GP) equation, which 
in this context is better known as the discrete nonlinear 
Schrodinger (DNLS) equation. By contrast, the trun- 
cated Wigner phase-space method evolves an ensemble 
of points according to the DNLS equation, and thereby 
takes into account the non-linear squeezing or stretching 
of the distribution. 

In the semiclassical treatment the quantum state in 
any moment is regarded as a "mixture" of wavefunctions 
ipi rather than a single V'i- It is worth noting that the 
stationary solutions of the DNLS equation are simply the 
fixed points of the Hamiltonian fiow. The small oscilla- 
tions obtained by linearization around these fixed points 
are the so-called Bogoliubov excitations. The typical 
oscillation frequency of the Bloch vector generally ap- 
proaches the classical frequency as N is increased keep- 
ing u fixed. However, in the vicinity of the separatrix 
convergence to the (vanishing) classical frequency is log- 
arithmically slow, as found via WKB quantization. 

Based on the ratio between the width of the semiclas- 
sical distribution for SCS and the width of the separatrix 
phase-space region, we find that the long-time Fringe Visi- 
bility of an initially coherent state in the Josephson inter- 
action regime, has a u/N dependent value (Fig. [7]). The 
functional dependence on u/N varies according to the 
preparation. In particular, whereas a Zero relative-phase 
preparation remains roughly Gaussian (i.e. phase-locked) 
throughout its motion, thereby justifying the use of MET 
for the description of Josephson oscillations around it, a 
Pi relative-phase SCS squeezes rapidly and its relative- 
phase information is lost [l^jH^- contrast, starting 
from fully separated modes, the phase distribution in the 
Josephson regime assumes a non-uniform profile, peaked 
at = 0, yielding a universal Fringe Visibility value of 

1/3 m. 

Focusing on two types of coherent preparations in 
the vicinity of the separatrix we find significant differ- 
ences in their M dependence on {u;N). The Pi SCS 
preparation (with vanishing population imbalance and 
a TT relative-phase) exhibits u dependent fluctuations, 
whereas the Edge SCS (having a comparable population 
imbalance but located elsewhere along the separatrix) ex- 
hibits N dependent fluctuations. Only in the latter case 
is the classical limit approached easily by taking large N 
at fixed u. 



To conclude, we have applied a semiclassical phase- 
space picture to the analysis of the one-particle coher- 
ence loss and buildup in the bosonic Josephson junction, 
described by the two-site BHH. The simplicity of the clas- 
sical phase space of the dimer allows for its semi-analytic 
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Appendix A: An L site system with N bosons 

A bosonic L site system has formally the same Hilbert 
space as that of L coupled (harmonic) oscillators, with 
an additional constant of motion TV that counts the to- 
tal number of quanta (rather than the total energy). 
The creation operator a| corresponds to a raising op- 
erator and the occupation of the ith site rii = a|ai 
corresponds to the number of quanta stored in the 
ith mode. The one-particle states of an L site sys- 
tem form an L dimensional Hilbert space. The set 
of unitary transformations (i.e. generalized rotations) 
within this space is the SU{L) group. For N parti- 
cles in L sites, the dimensionality of the Hilbert space 
is dim(Af) = {N+iy./[{L-iy.{N-L+2)\]. For example 
ioT L — 2 we have dim(iV) = A^-|-l basis states \ni,n2) 
with ni + n2 — N. We can rotate the entire system using 
dim(iV) matrices. Thus we obtain a dim(iV) representa- 
tion of the SU{L) group. By definition these rotations 
can be expressed as a linear combination of the SU{L) 
generators J^, and they all commute with the conserved 
particle number operator A'". The many-body Hamilto- 
nian % may contain "nonlinear" terms such as J'^^ that 
correspond to interactions between the particles. Accord- 
ingly, for an interacting system, % is not merely a rota- 
tion. However, % still commutes with N ^ maintaining 
the fixed-number dim(A'^) subspace. 

In the semiclassical framework the dynamics in phase 
space is generated by the Hamilton equations of motion 
for the action-angle variables hi and tpi. These are the 
"polar coordinates" that describe each of the oscillators. 
It is common to define the complex coordinates 



(Al) 



(representing a single point in phase space). This is the 
classical version of the destruction operator a^. The equa- 
tion for \E'i is the discrete nonlinear Schrodinger (DNLS) 
equation: 



dt J - . 2 

which is the space-discretized version of 



(A2) 



i ; 

dt 



Vix) 



-g,\^{x)\'^^\/'U{x) (A3) 
zm J 



This looks like the Gross-Pitaevskii (GP) equation, but 
strictly speaking it is not the GP equation. The GP equa- 
tion is the outcome of a mean-field theory (MFT): it is 



not an equation for ^'i, but an approximated equation 
for the mean-field 5*^. We further clarify this point in 
the next paragraph. 

Within the framework of the semiclassical treatment 
the quantum state is described as a distribution of points 
in phase space. This approach goes beyond the con- 
ventional MFT approximation: MFT essentially assumes 
that the state of the system remains coherent throughout 
its evolution. Such a state corresponds to a Gaussian-like 
distribution is phase space ( "minimal wave-packet" ) and 
accordingly it is fully characterized by "center of mass" 
coordinates {Jp^,ni) that are defined through the mean 
field 



Mean Field 



(A4) 



(representing the center of a wave-packet). To the extent 
that the MFT assumption can be trusted, the equation 
of motion for is the DNLS / GP. Indeed, if u = 
there is no nonlinear spreading and the MFT descrip- 
tion becomes exact. The approximation holds well in 
the weak-interaction Rabi regime and for some regions 
of phase-space (e.g. around the ground state) in the 
strong-interaction Josephson regime, but in the nonlin- 
ear regions of phase-space MFT becomes too crude to 
provide an accurate description of the dynamics. 



Appendix B: The structure of phase space 

The energy contours in a representative case are dis- 
played in Fig. [H Considering a section along the big 
circle = 0,7r, it is convenient to take — tt <9 < -|-7r 
instead of < < tt. Along this section the energy is 

E{e) = {NK/2)f{e) where 



I{0) 
fid) 

no) 



e cos 9 — sin 9 



(Bl) 



-iisin(20) +£sin6'- cos6', (B2) 



= — ^tcos(2^^) + £cos0 + sin( 



(B3) 



All extremal points of the Hamiltonian are located 
along this section and are determined by the equation 
f'{9) = 0. The number of solutions depends on u and on 
the bias e. 

For u < 1, there are two fixed points: a minimum at 0_ 
at the bottom of the sea and a maximum at the opposite 
point 9+. For e = the two fixed points are: 



7r/2 
-7r/2 



[i^=0 point], 
['f—TT point], 



(B4) 
(B5) 



with corresponding energies E± — ±{N/2)K. 

For li > 1, provided |e| < ec, the upper fixed point 
bifurcates into a saddle point 9y^ and two stable maxima 
6*1 and 02- The value of the critical bias Ec is derived in 
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the next paragraph. For e = the four fixed points are: 



9- = +7r/2, 
0x = -n/2, 
9i_2 = — arcsin(l/u), 



(B6) 
(B7) 
(B8) 



with corresponding energies E^, and E+ that 

are given by Eqs. dlOIlT^ . We can also determine 
the borders of the separatrix by solving the equation 
/(^i',2') = /(^x)- Thus we get the following expressions 
for the outer borders of the islands 



arcsin(l — (2/u)). 



(B9) 



The value of the critical bias Eq. © is obtained by 
solving the equation f'{6) — together with f"{6) — 0. 
An equivalent set of questions is 

sin(6')/"(6') - cos{9)f'{9) = 1 + usm^{9) = 0, (BIO) 
cos{9)f"{9) + sin{9)f'{9) =e-u cos^ 9 = 0. (Bll) 

leading to the solution 



= - arcsin((l/u)^/^ 

3/2 



(B12) 

(^2/3 _ (B13) 



For Ec = £, one island has zero area while the other has a 
critical area Ac. For completeness we derive an estimate 
for Ac, that we have used in [l^. The island is defined 
by the equation f(9, ip) = const., where 

f{9,ip) = ^u{cos9)'^ - ecos9 - sm9cos(p. (B14) 

The critical island is defined by the equation f{9, ip) = fc 
where according to the previous paragraph 



fc 



f{9.\ec 



2 2 



(B15) 



The outer turning point 9y is identified as the second 
root of the equation f{9) = fc- Defining s = cos6' the 
equation takes the form 

_ [ecu]s^ + [1 - fcu+el]s^ + [2ecfc]s + [fl-l] = 0. 
One root with double degeneracy is obviously 

= COS(0,) = (l-M"2/3)l/2^ (B^g) 

Then we can get the third root by solving a quadratic 
equation, leading to 



(B17) 



The area of the critical island can be found numerically, 
and a very good approximation is s^, namely 



47r (1-^-2/3)3/2. 



(B18) 



Appendix C: The Wigner function of a spin 

The Stratonovich- Wigner- Weyl correspondence 
(SWWC) associates with any Hermitian operator of a 
spin^^', a real sphere'^-'-' function A„(ri). Both spaces 
have the same dimension {2j+l)'^, and the association 
is one-to-one. This appendix provides the reasoning and 
the practical formulas following [4^ 13] ■ 

The inner product of two operators is defined as 
trace[A'l'B]. An orthonormal-like set of projector-like op- 
erators, that are knows and the Stratonovich- Weyl oper- 
ators, can be defined, such that 



trace P^P^ 



i, (CI) 

Sj{n-n') (C2) 



where ft = {9,(p), and the "delta" function is 



2j e 



5j{n-n') = J2 J2 Y''"'{n)Y'^"'*{n') (ca) 



£=0 m=-£ 



Accordingly any operator A can be represented by the 
phase-space function 



A„{n) = trace[P" p]. 



(C4) 



From the above definition it follows that the inner prod- 
uct of two Hermitian operators is given by a phase-space 
integral over the product of the corresponding Wigner 
functions as in Eq. ([32]) . 

The actual construction of the P^ is not a simple task. 
One procedure is to start with the non-orthonormal set of 
coherent state projectors |0)(0| and to perform an "or- 
thogonalization" , very similar to that employed in con- 
dense matter physics for the purpose of defining a Wan- 
nier basis. The final result is conveniently expressed as 



2j 



- 2j 



J2 Y^"'{n) f''". 



(C5) 



which is an orthogonal transformation over the non- 
Hcrmitian mulitpole operators T'™ defined as 



f""= y (-l)^-™'x/2]TTf ^ , ^ )|m')(m" 
^ ' \ —TO ni TO ' ' ' ^ 



We use here the Wigner 3j symbols. Conse- 
quently we can represent any operator A either by 
Aim = trace[(f '™)ti] or by A„{Vl) = trace[P"i]. In 
particular it follows that p = |'0)('0I is represented by 



l,m 



where 



m' .7n" 



2j+l 



j I j 
-to' to to'' 



(C6) 
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For the coherent state p = |i, j | the only term in the This function resembles a minimal Gaussian as one can 
sum is m' — m" = j, and therefore all the Cim are zero see from the example in Fig. [21 
except m — 0. Consequently, the Wigner function is 

I 

Appendix D: LDOS semiclassical calculation 



Using the definition Eq. (pS)) we can calculate the LDOS in the semiclassical approximation using Eq. ([37]). The 
integral that should be calculated in the case of Zero {cp- ~ 0) and Pi {ip^ — tt) preparations is 



P(£;) = uj{E) 



dipdn 
2tt 



NK 



■ coa{(p) — E 



1 



■ cxp 



2a2 



2&2 



(Dl) 



with E — Eij. Most of the contribution comes from the vicinity of the fixed point, so we can use a quadratic 
approximation: 



V{E) uj{E) 



dipdn 
2-K 



1 



(5 C/n^ ± -NKlp^ - (E-E^^^) — exp 



1 



ab 



'2a2 



n 

262 



Below we derive the following results for the Zero and Pi preparations: 



P(£') = 2 exp 



f-L 



4 



{E-E^) 



lo 



4 



1 

ivcJ 



p(^;) . i f ^ 



exp 



1-^^ iE- 



K NU J 



(E-E^ 



(E-E^) 



±^^\lE- 



K NU 



E-E^ 



(D2) 

(D3) 
(D4) 



In order to simplify notations we define bellow e as the difference E—E^ or E—E^. 

In order to estimate the integral in the Zero case we change to polar coordinates r and — tt < t < it: 

LP = [Ae/NKY''^ rcos(t), 
n = [e/UY'"^ rsin(t), 



leading to 



p(£;) 



Wo 

1 

2^ 
1 

TT 

exp 



rdrdt 2 



27r loq 

dt — exp 
ah 



■ cxp 



-{Ae/NK) 



+7r 



dt exp 



1 



K 2NU 



2e e 
K 2NU 

2Io 



2e e 
K ^ 2NU 
1 



2a2 
262 

cos(2t) 



262 



K 2NU 



(D5) 
(D6) 



(D7) 
(D8) 
(D9) 
(DIO) 



It is important to realize that the semiclassical evaluation in the Zero case is valid only if the contour lines of the 
Gaussian intersect the contour lines of the eigenstates transversely. Else, the Airy structure of the eigenstates should 
be taken into account, or perturbation theory rather than semiclassics should be used. The asymptotic behavior of the 
Bessel function is Io(a;) « exp(x)/V27rx and hence P{E) ^ exp[~e/{NU)]. One observes that the tails of the LDOS 
reflect the relatively slow decay of the Gaussian tails in the n direction. 

In order to estimate the integral in the Pi case we change to the coordinates r and — oo < t < oo: 



ip = [Ae/NK]^/^ rsinh(f) 
n= [e/C/]i/2 rcosh(t) 



or cosh(t) for e < 0, 
or sinh(t) for e < 0, 



(Dll) 
(D12) 
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leading to 



F{E) 



rdrdt 2 



27r ujq 



1 fujiE) 

IT 



1) 



1 



ab 



exp 



(sinh(t))= 
2^2 



(|,|/f/)l£^!i£))!^2 



dt — exp 
ab 

^0 ) 



-(4kl/ivx) - (H/t/) 



dt exp 



± 



1 



Wo 



exp 



1 



X 2NU 



2N 
Ko 



2N 
if 



2iVL/ 



262 

(C0sh(t))2 
2fe2 

cosh(2i) 



2 1 

X ^ 27VC7 



(D13) 
(D14) 

(D15) 

(D16) 



Note that without the Gaussian the resuh of the integral should be unity reflecting the proper normalization of 
the microcanonical state. The upper cutoff of the t integration reflects the finite size of phase space. The cutoff 
prevents the singularity in the limit e — >■ 0. The Bessel function expression is obviously valid only outside of the 
cutoff-affected peak. The asymptotic behavior of the Bessel function is Ko{a;) ~ exp(— a::)/-\/27rx, and hence we get 
P(£;) ~ exp[-|e|/(iVC/)] for e > and V{E) - exp[-|e|/(if/4)] for e < 0. One observes that the tails of the LDOS for 
e < reflect the relatively rapid decay of the Gaussian tails in the direction. 

Finally we point out that the calculation of the LDOS for the Edge and for the TwinFock preparations is much 
easier. In the former case the delta function of p^^^{n^ ip) merely replaces the n coordinate in p^^'>{n, ip) by a constant 
proportional to E—E^, leading to a Gaussian for P{E). Similarly for the TwinFock preparation: 



PiE) 



u{E) 



'^'P^'' 5(un'~^co.{^)-E] 5{n) oc 



u(E) 



27r " V" 2 J " {NKl2)\sin{if)\_ 

The later expression should be calculated for {NK/2) cos,{tf) ~ E leading to Eq. (PT|) in the main text. 

I 



(D17) 



Appendix E: Rabi-Josephson oscillations, MFT and 
semiclassical perspectives 

MFT.— For u = 0, the Hamiltonian ^ merely gen- 
erates rotations in phase space. Consequently, coherent 
states remain Gaussian-like throughout their motion, re- 
sulting in Rabi oscillations of the population between the 
sites. This is the case where MFT gives exact results. 
For u 7^ 0, MFT maintains the assumption that the sys- 
tem remains in a coherent state at any instant during its 
evolution and that propagation only serves to displace 
this Gaussian distribution on the Bloch sphere. Thus, 
MFT corresponds to the classical dynamics of a "point" 
in phase space. In particular, this implies that the occu- 
pation statistics is binomial at any time. 

Nonlinear effects cannot be neglected once the wave- 
packet is stretched along the phase-space energy contour 
lines. Then MFT no longer applies, but the semiclassi- 
cal approximation still works well for the description of 
such squeezing. For example, by projecting the evolving 
semiclassical distribution onto n we obtain the number 
distribution Pt (n) ^25j] whose line shape reflects the phase 
space structure, with caustics at the borders of forbidden 
regions (see Fig. [3]). 

Semiclassics. — For m = 0, all trajectories have the 
same topology and frequency and a Gaussian-like wave- 
packet that is launched (say) at the NorthPole executes 
Rabi oscillation between the two wells. If u is non-zero 



but small, then all the trajectories still have the same 
topology, but uj{E^) is v dependent. The anharmonic be- 
havior is due to the spectral stretch Awoa^- Accordingly, 
we distinguish between the harmonic stage (t < l/Awosc) 
and the anharmonic stage {t > I/Awobc) in the time evo- 
lution. 

For u > 1, which we call "the Josephson regime", a 
separatrix emerges and accordingly there are two types 
of stable coherent oscillations, depending on the initial 
population imbalance iJJj] : Small oscillations around 
the ground state, lying in the bottom of the sea; and 
self-trapped oscillations around the top of either island. 

For u > N"^ which we call "the Fock regime" , the area 
of the sea becomes less than a Planck cell, and therefore 
effectively disappears. Our main interest in this work is 
the unstable motion along the separatrix for 1 <C u <C iV^ 
[2^ - 12^ . Such motion emerges for the Pi and for the Edge 
preparations of Fig. [2j 

Beyond. — Quantum effects that are ignored by the 
semiclassical approximation are anharmonic beats, long- 
time recurrences, and the possibility of tunneling between 
the two islands. Associated with the recurrences are the 
long-time fluctuations, both in the occupation and in the 
Fringe Visibility, as discussed in Section VIII. 



14 



The [arXiv: 1001.2120 version of this paper includes 3 ex- 
tra pedagogically oriented appendices: Bosons in A'^ site [27 
system and MFT; Detailed analysis of phase space geom- 
etry in the A'^=2 case; The definition of a Wigner function 
of a spin; Concise discussion of the Rabi-Josephson oscil- [28 
lations in the semiclassical perspective. 

D. Jaksch C. Bruder, J. I. Cirac, C. W. Gardiner, and P. 
ZoUer, Phys. Rev. Lett. 81, 3108 (1998). [29 
M. Greiner, O. Mandel, T. Esslinger, T. W. Hansch, and 
I. Bloch, Nature (London) 415, 39 (2002). [30 
Y. Makhlin, G. Schon, and A. Shnirman, Rev. Mod. 
Phys. 73, 357 (2001); R. Gati and M. K. Oberthaler, [31 
J. Phys. B 40, R61 (2007). [32 
Gh-S. Paraoanu et al, J. Phys. B: At. Mol. Opt. Phys. 

34, 4689 (2001); A. J. Leggett, Rev. Mod. Phys. 73, 307 [33 
(2001). 

J. Javanainen, Phys. Rev. Lett. 57, 3164 (1986). 
F. Dalfovo, L. Pitaevskii, and S. Stringari, Phys. Rev. A [34 
54, 4213 (1996). 

I. Zapata, F. Sols, and A. J. Leggett, Phys. Rev. A 57, [35 
1050 (1998). 

A. Smerzi, S. Fantoni, S. Giovanazzi, and R. S. Shenoy, [36 
Phys. Rev. Lett. 79, 4950 (1997). 

F. S. Cataliotti, S. Burger, C. Fort, P. Maddaloni, F. [37 
Minardi, A. Trombettoni, A. Smerzi, and M. Inguscio, [38 
Science 293, 843 (2001). 

M. Albiez, R. Gati, J. Foiling, S. Hunsmann, M. Oris- [39 
tiani, and M. K. Oberthaler, Phys. Rev. Lett. 95, 010402 
(2005). [40 
S. Giovanazzi, A. Smerzi, and S. Fantoni, Phys. Rev. [41 
Lett. 84, 4521 (2000). 

S. Levy, E. Lahoud, I. Shomroni, and J. Steinhauer, Na- [42 
ture 449, 579 (2007); C. A. Sackett, Nature 449, 546 
(2007). 

A. J. Leggett and F. Sols, Found. Phys. 21, 353 (1998). [43 

E. M. Wright, D. F. Walls and J. C. Garrison Phys. Rev. 
Lett. 77, 2158 (1996). 
J. Javanainen and M. Wilkens, Phys. Rev. Lett. 78, 4675 [44 
(1997); Phys. Rev. Lett. 81, 1345 (1998). 

M. Greiner, O. Mandel, T. W. Hansch, and I. Bloch, 
Nature 419, 51 (2002). [45 

G. -B. Jo et al, Phys. Rev. Lett. 98, 030407 (2007). 
T. Schumm et al, Nat. Phys. 1, 57 (2005). [46 
S. Hofferberth et al. Nature (London) 449, 324 (2007). 

A. Widera, S. Trotzky, P. Cheinet, S. FoUing, F. Gerbier, [47 
I. Bloch, V. Gritsev, M. D. Lukin, and E. Demler, Phys. 
Rev. Lett. 100, 140401 (2008). [48 
A. Vardi and J. R. Anglin, Phys. Rev. Lett. 86, 568 [49 
(2001); J. R. Anglin and A. Vardi, Phys. Rev. A 64, 
013605 (2001). [50 
Y. Khodorkovsky, G. Kurizki, and A. Vardi, Phys. Rev. [51 
Lett. 100, 220403 (2008); Phys. Rev. A 80, 023609 [52 
(2009). 

E. Boukobza, M. Chuchem, D. Cohen, and A. Vardi, [53 
Phys. Rev. Lett. 102, 180403 (2009). 

K. Smith-Mannschott, M. Chuchem, M. Hiller, T. Kot- [54 
tos, and D. Cohen, Phys. Rev. Lett. 102, 230401 (2009). 
C. W. Gardiner and P. ZoUer, Quantum Noise, 2nd ed. 



(Springer, New York, 2000). 

M. J. Steel, M. K. Olsen, L. I. Phmak, P. D. Drummond, 
S. M. Tan, M. J. CoUett, D. F. Walls, and R. Graham, 
Phys. Rev. A 58, 4824 (1998). 

A. Sinatra, C. Lobo, and Y. Castin, Phys. Rev. Lett. 87, 
210404 (2001); J. Phys. B: At. Mol. and Opt. Phys. 35, 
3599 (2002). 

I. Carusotto, Y. Castin, and J. Dalibard, Phys. Rev. A 
63, 023606 (2001). 

S. E. Hoffmann, J. F. Corney, and P. D. Drummond, 

Phys. Rev. A 78, 013622 (2008). 

A. Polkovikov, Phys. Rev. A 68, 053604 (2003). 

L. I. Plimak, M. Fleischhauer, M. K. Olsen, and M. J. 

Collett, Phys. Rev. A 67, 013812 (2003). 

P. Deuar and P. D. Drummond, J. Phys. A: Math and 

Gen. 39, 1163 (2006); J. Phys. A: Math and Gen. 39, 

2723 (2006); 

L. Isella and J. Ruostekoski, Phys. Rev. A 74, 063625 

(2006) . 

M. Hiller, T. Kottos, and T.Geisel , Phys. Rev. A 73, 
061604(R) (2006); Phys. Rev. A 79, 023621 (2009). 
S. L. W. Midggley, S. Wiister, M. K. Olsen, M. J. Davis, 
and K. V. Kheruntsyan, Phys. Rev. A 79, 053632 (2009). 
P. Deuar, Phys. Rev. Lett 103, 130402 (2009). 

F. Trimborn, D. Witthaut, and H. J. Korsch, Phys. Rev. 
A 77, 043631 (2008); Phys. Rev. A 79, 013608 (2009). 
K.W. Mahmud, H. Perry, and W.P. Reinhardt, Phys. 
Rev. A 71, 023615 (2005). 

Franzosi et al. Int. J. Mod. Phys. B 14, 943 (2000). 

E. Boukobza, D. Cohen, and A. Vardi, Phys. Rev. A 80, 

053619 (2009). 

E.M. Graefe, H.J. Korsch, Phys. Rev. A 76, 032116 

(2007) ; D. Witthaut, E. M. Graefe, and H. J. Korsch, 
Phys. Rev. A 73, 063609 (2006). 

G. S. Agarwal, Phys. Rev. A. 24, 2889 (1981). J. P. Dowl- 
ing, G. S. Agarwal,W. P Schleich, Phys. Rev. A 49, 4101 
(1994). 

J.C. Varilly and J.M. Gracia-Bondia, Annals of Physics 
190, 107 (1989). C. Brif and A. Mann, J. Phys. A 31, 
L9 (1998). 

A. Polkovnikov, S. Sachdev, and S.M. Girvin, Phys. Rev. 
A 66, 053607 (2002). 

E. Altman and A. Auerbach, Phys. Rev. Lett. 89, 250404 
(2002). 

A. K. Tuchman, C. Orzel, A. Polkovnikov, and M. Kase- 
vich, Phys. Rev. A 74, 051601 (2006). 

B. Wu and Q. Niu, Phys. Rev. A 61, 023402 (2000). 

J. C. Eilbeck, P. S. Lomdahl, and A. C. Scott, Physica 
D 16, 318 (1985). 

R. L. Stratonovich, Sov. Phys. JETP 31, 1012 (1956). 

D. Cohen and T. Kottos, Phys. Rev. E 63, 36203 (2001). 

E. Boukobza, M. G. Moore, D. Cohen, and A. Vardi, 
Phys. Rev. Lett. 104, 240402 (2010). 

M. Rigol, V. Dunjko, and M. Olshanii, Nature 452, 854 

(2008) . 

Quantum Optics by M.O. Scully and M.S. Zubairy, 
(Cambridge University Press, 1997), p. 201. 



15 




FIG. 1: (Color online) Contour lines for u > 2. Sea levels are colored blue, Island levels are colored green, and the Separatrix is 
colored red (left panel). Energy spectrum for A'' = 20 and u = 10. WKB energies (red x) are compared with exact eigenvalues 
(blue +). Dashed lines indicate slopes ojj for low energies, ojx for near-separatrix energies, and u!+ for high energies (right 
panel) . 



1 




FIG. 2: (Color online) An illustration of TwinFock (n = 0) preparation (left), and of Pi ("tt"). Zero ("0") and Edge ("e") 
preparations (right) using Wigner plots on a sphere. The left panel is a 3D plot, while the right panel is a Mercator projection 
of the sphere using [tp, n) coordinates. 
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FIG. 3: (Color online) The evolving quantum state of A'^ = 40 bosons with u — 5 for TwinFock (n = 0) preparation. Here and 
below the units are such that K — 1. The time is t — 4. (a) Wigner function of the evolved quantum state, (b) Semiclassical 
evolved state, (c) Occupation statistics, with the semiclassical result shown as dashed-dotted line. 
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FIG. 4: (Color online) The variation of Sx{t) with time for A*'=40 particles with u — 5, for Zero (a), Pi (b), Edge (c), and 
TwinFock (d) preparations. Note the different vertical scale. The dashed-dotted lines are based on semiclassical simulation. 
Note that the fluctuations of a semiclassical preparation always die after a transient, which should be contrasted with both the 
classical (single trajectory) and quantum (superposition of Af > 1 eigenstates) behavior. 
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FIG. 5: (Color online) The LDOS (left) and the spectral content of the fluctuations (right) of TV = 500 bosons with m = 4, for 
Zero, Pi, Edge, and TwinFock preparations (top to bottom). The horizontal axes are E — and ui/uj. The lines in the LDOS 
figures are based on a semiclassical analysis (see text), while the circles are from the exact quantum calculation. Note that due 
to the mirror symmetry of the Zero preparation the expected frequency should approach 2a; j, while for the Pi preparation it 
is bounded from below by 2u}^ (both frequencies are indicted by vertical dashed lines). Note also the outstanding difference 
between the spectral support of Zero and Pi preparations compared with continuous-like support in the case of Edge and Fock 
preparations. 
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FIG. 6: (Color online) The participation number M as determined from the LDOS for N = 100 (o), 500 (□), and 1000 (o) 
particles. The left panel contains the Zero (lower set in blue) and Pi (upper set in red) preparations, while the Edge preparation 
is presented in the right panel. Note the different vertical scale. In the crudest approximation we expect in the Edge case 
M ~ iV^/^, while in the Pi case M < N'^^^ as long as (u/N) < 1 
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FIG. 7: (Color onhne) The mean frequency of the Sx{t) oscillations versus u/N for for A'' = 100 (o), 500 (□), and 1000 (<>) 
particles. The preparations are (upper to to lower sets of data points): Zero (blue), Edge (magenta), and Pi (red). The 
doubled Josephson frequency 2ujj is marked by a dashed blue line. The theoretical predictions of Eq. (|42[l . doubled due 
to mirror symmetry, are represented by red and magenta dashed lines, while Eq. (|43() for u/N S> 1 is represented by black 
dash-double-dotted line (there is one fitting parameter as explained there). 
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FIG. 8: (Color online) Left: The long-time average of S^t) versus u/N for = 100 (o), 500 (□), and 1000 (o) particles. The 
preparations are (upper to lower sets of data points): Zero (blue), TwinFock (black). Edge (magenta), and Pi (red). The 
symbols are used for the quantum results and the dashed lines are the semiclassical prediction for fourty particles. Note that 
the scaling holds only in the Josephson regime 1 <^ u <^ N^ , and therefore, for a given u/N range, becomes better for large A*'. 
Right: The long time RMS of Sx{t) for the three coherent preparations (lower to upper sets): Zero (blue), Edge (magenta), and 
Pi (red). The implied A''^''" scaling based on Eq. ^) is confirmed. In the inset, the RMS of Sx{t) for Edge (A) and Pi (v) 
preparations is plotted versus A'' while u = 4 is fixed. The dashed lines are power-law fits that nicely agree with the predictions 
of Eq. GOl). 



